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Foreword 


This  work  was  performed  during  the  period  August  1967  through 
December  1969  under  U.S.  Army  Natick  Laboratories  Contract  No. 
DAAG-17-67-C-0189  for  the  Department  of  the  Army  Project  No. 
1M121401  D195  entitled  "Exploratory  Development  of  Airdrop  Systems" 
Task  13  -  Impact  Phenomena,  The  program  is  a  part  of  continuing 
investigation  directed  toward  obtaining  improved  energy  dissipater 
materials  for  airdrop  landing  shock  mitigation  and  a  better  under¬ 
standing  of  the  response  of  airdroppable  materiel  to  airdrop  im¬ 
pact  phenomena. 
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ABSTRACT 


The  accelerations  and  displacements  in  a  complex  structure 
subjected  to  an  impact  loading  are  computed  by  treating  the 
structure  as  a  lumped  parameter  system.  A  mathematical  model  of 
the  system  consists  basically  of  discrete  masses  linked  by  weight' 
less,  elastic  beams  with  the  appropriate  stiffnesses,  areas,  and 
moment-of-inertia  properties.  By  specifying  a  proper  set  of  in¬ 
dependent  coordinates  through  which  the  motion  of  these  lumped 
masses  are  uniquely  described,  and  by  writing  equations  of  motion 
in  terms  of  these  coordinates,  a  set  of  equations  is  derived 
which  represents  the  motion  of  any  part  of  the  model  during  im¬ 
pact.  Using  the  Runge-Kutta  numerical  method,  and  a  digital  com¬ 
puter,  these  equations  are  solved. 

A  physical  model  of  the  lumped  parameter  system  was  built 
and  cushioned  with  paper  honeycomb.  Displacements  and  accelera¬ 
tions  at  some  points  in  this  model  were  measured  and  compared 
with  computed  results.  Agreement  is  satisfactory. 


IMPACT  ON  COMPLEX  STRUCTURES 


1.  Introduction 

This  study  has  two  objectives;  (1)  to  learn  how  to 
realistically  model  a  vehicle  using  lumped  masses  and 
springs  and  to  relate  such  computed  quantities  as  dis¬ 
placement  and  accelerations  to  corresponding  quantities 
in  the  prototype;  and  (2)  to  prepare  a  computer  code  for 
doing  the  necessary  computations,  and  which  is  sufficiently 
flexible  to  serve  as  a  design  tool. 

The  report  consists  of  eight  sections  as  foDlows: 

Lumped  mass  model 
Cushioning  system  design 
Analytical  study 
Numerical  solution 
Digital  computer  code 
Experimental  work 
Discussion 
Conclusion 

2.  Lumped  Mass  Model 

A  complex  structural  body  can  be  represented  by  a  system  of 
discrete  lumped  masses  connected  by  massless  beams  with  appropri¬ 
ate  flexural  stiffnesses.  In  airdrop  practice  the  bending  mo¬ 
ments  and  stresses  in  the  structure  are  minimized  by  dividing  the 
structure  conceptually  into  free  bodies,  cushioning  each  part  in¬ 
dependently  to  provide  zero  rigid  body  rotation  and  equal  rigid 
body  translatory  acceleration  for  each  part.  To  accomplish  this, 
the  cushioning  forces  must  be  distributed  over  the  structure  in 
proportion  to  the  weight. 

The  model  shown  in  Figs .  1  and  2  may  be  regarded  as  a  highly 
simplified  representation  of  the  M37  truck.  The  mass  distribution 
of  the  truck  is  roughly  approximated,  and  the  various  elastic  beam 
elements  simulating  the  structural  members  are  assumed  to  have 
equal  stiffness  properties.  No  damping  has  been  considered.  Mass 
1  can  be  regarded  as  the  motor  and  mass  2  the  load  in  the  vehicle. 
Masses  3  and  4  hanging  on  springs  below  the  main  masses  are  like 
the  wheels.  Two  small  masses,  5  and  6,  are  attached  in  cantilever 
fashion  to  the  main  masses.  These  two  do  not  represent  anything 
in  particular.  They  are  included  to  provide  information  on  how 
much  such  masses  influence  the  motion  of  the  mam  mass,  and  how 
the  accelerations  of  the  small  masses  are  modified  by  their  loca¬ 
tions  within  the  structure.  Mass  7  (Fig.  2),  which  is -3$  of  the 


1 


■X 


total  model  mass,  has  also  been  included  to  show  how  a  mall 
mass  not  attached  directly  to  another  mass  may  influence  the 
motion  of  the  main  masses. 

The  mathematical  model  that  corresponds  to  the  physical 
model  is  shown  in  Fig.  3.  If  mass  7  is  connected  directly 
to  the  beam  its  very  small  mass  in  comparison  with  the 
other  masses  results  in  a  large  acceleration  (a  =  F/m)  which 
causes  computational  difficulties.  To  avoid  these  diffi¬ 
culties^  a  small  beam  element  is  used  to  connect  this  mass 
to  the  large  beam. 

3.  Cushioning  System  Design 

a.  Forces 

The  overall  magnitudes  of  the  cushioning  forces  applied 
to  the  model  as  shown  in  Fig.  3  are  determined  from  the 
acceleration  level  desired  during  impact.  For  a  stated 
acceleration  level  of  G,  measured  in  units  of  acceleration 
of  gravity,  the  total  force  applied  is 

EF  =  FA  +  Ffi  +  Fc  +  FD  +  FE  *  (U1  +  W2  +  W3  +  W4)  (G  +  l) 

To  reduce  bending  moments  in  the  long  structural 
spans  between  masses,  a  vehicle  that  is  to  be  cushioned 
is  subdivided  into  free  bodies,  and  force  and  moment 
balances  are  carried  out  for  each  section.  Each  free  body 
has  to  be  a  statically  determinate  structure  and  no  force 
~r  moment  is  transferred  to  this  free  body  through  the 
point  of  cut.  If  a  cushioned  point  belongs  to  several 
free  bodies  then  the  cushioning  force  at  that  point  will 
be  the  sum  of  cushioning  forces  of  all  free  bodies  at  that 
point.  For  this  analysis,  the  model  has  been  divided  into 
three  sections.  Free  bodies  of  these  three  sections  are 
shown  in  Fig.  *1.  The  cushioning  in  the  idealized  mathe¬ 
matical  model  provides  the  forces  shown. 

b.  Cushioning  Material  Characteristics 

The  material  composing  the  cushioning  system  is 
idealized  by  assuming  that  it  crushes  at  a  constant  strength 
and  that  it  has  no  resilience.  This  is  an  idealization  of 
the  characteristics  of  paper  honeycomb.  The  computer  pro¬ 
gram  at  the  present  time  includes  only  a  constant  crushing 
force.  Time  varying  forces  will  be  considered  later.  To 
obtain  data  applicable  to  the  cushions  used  in  the  labora¬ 
tory  model  nine  3  inch  x  3  inch  x  3  inch  paper  honeycomb 
pads  were  tested  with  a  drop  height  of  38  inches  and  a  mass 


VEHICLE  SUBDIVISION 

SECTION  (A) 


SECTION  (B) 


F,*  Fc*-^Wi(6+l) 


SECTION  (C) 


Fe-  W4(G+I) 


FIG.  4  FREE  BODY  DIAGRAMS 


weight  of  220  pounds.  A  typical  str®®s7®tJ^n  ib/ft2* 
shown  in  Pig.  5.  The  average  strength  is  kOOO  lb/ft  . 
This  value,  due  to  the  small  size  of  the  pad,  is  much 
smaller  than  the  normal  value  for  paper  honeycomb  (6000 
psf).  Based  on  the  assumed  constant  crushing  force  oi 
the  paper  honeycomb  the  design  of  the  cushioning  system 
is  as  follows: 


FIG.  5  TYPICAL  STRESS  -  STRAIN  RECORD  FOR  3  in. 
X  3  in.  X  3  in.  PAPER  HONEYCOMB  SAMPLES 


c.  Design  Procedure 

First  select  a  design  acceleration  leyel  G  and  use 
Newton’s  law  to  compute  cushioning  force.  • 


where  Fp  ■  the  resultant  force  applied  to  the  mass 
M  ■  —  ■  Mass 
W  ■  weight  of  mass 

a  *  acceleration  of  the  center  of  gravity  of  the  mass 
g  3  32.2  ft/sec2 

Equation  (1)  can  be  rewritten  as 

Fr  »  WG - (2) 

Q 

where  G  *  —  *  a  dimensionless  number  which  indicates  the  ac 
K  celeration  in  "g's." 

The  relation  between  effective  cushioning  area  A  and  the  as 
sumed  acceleration  level  G,  and  the  crushing  strength  S 
is  derived  as  follows,  referring  to  Fig.  6. 


FIG.  6  IDEALIZED  CUSHIONING  SYSTEM 

F  -  W  »  W(|)  -  WG 

O 

or 

F  3  SA  3  W(G  +  1) 

and 

A  3  g(G  +  1) - (3) 

Computation  of  the  required  thickness  z  of  the  cushion 
is  as  follows. 


8 


where  V  «  required  volume 


The  volume  will  depend  on  the  amount  of  energy  to  be  dissi¬ 
pated,  which  will,  in  turn,  depend  on  the  impact  velocity 
or  the  drop  height  h. 

The  potential  energy  U  of  the  package  will  be 

U  =  W(h  +  ez) - (4) 

where 

c  *  design  strain. 


Since  the  stress-strain  curve  for  the  cushion  material 
is  assumed  to  be  approximated  by  a  rectangle,  the  amount 
of  energy  per  unit  volume  which  the  pad  will  absorb  is 


The  required  volume  •’s  therefore 
„  _  U  W(h  +  ez) 


and  the  necessary  thickness  is 

z  -  & - - - 

Oe 


The  area  and  thickness  of  the  required  cushion  stack 
thus  can  be  obtained  from  Equations  (3)  and  (6).  As 
stated  previously,  the  cushioning  forces  have  been  assumed 
constant  so  long  as  the  velocity  at  the  cushioned  points 
remains  negative,  indicating  downward,  hence  compressing 
movement.  However,  as  soon  as  the  velocity  becomes  zero 
or  changes  sign  this  cushion  force  will  vanish. 

4.  Analysis  of  the  Model 


a.  Equations  of  Motion 


The  general  form  of  the  equations  of  motion,  for  the 
model  in  matrix  notation,  is  . 


CM]  {<3)  +  [K]  {q}  -  {Q} - (7) 

where  [M]  represents  the  generalized  mass  matrix,  [K]^ 
represents  the  generalized  stiffness  matrix,  and  { Q } 
represents  the  generalized  force  matrix.  The  vector  {q} 


represents  the  independent  generalized  coordinate  system 
which  is  sufficient  to  uniquely  describe  the  positions  of, 
and  the  motion  of  the  entire  model,  while  is  the  vector 
describing  the  acceleration  of  the  model  along  these  coor¬ 
dinates.  The  generalized  coordinates  for  this  model  are 
shown  in  Pig.  3. 

b.  Element  Stiffness 

The  stiffness  matrix  for  the  individual  elements  is 
determined  from  elementary  beam  deflection  theory  and  the 
equations  representing  the  elastic  force-deflection  proper¬ 
ties  of  the  elements  are  combined  in  the  following  matrix 
expression. 

<F>  -  [K]  (u) 
where 


{F}  *  Force  vector 
[K]  *  Stiffness  matrix 
{u}  =  Displacement  vector 

The  coordinates  of  a  beam  element  and  a  spring 
element  are  shown  in  Fig.  7. 

The  stiffness  matrix  for  the  spring  element  is 


^K^spring  » 


k 


K 

o  — — VWVAA 


(a)  BEAM  ELEMENT 


(b)  SPRING  ELEMENT 


FIG.  7  ELEMENT  COORDINATES 
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The  element  stiffness  matrix  for  the  beam  is 


CR1  beam 


El 


12  -6L  -12  -6L 

-  6L  iJL2  6L  2L2 
-12  6L  12  6L 

-  6L  2L2  6L  *JL2 


The  generalized  stiffness  matrix  [K]  of  the  struc¬ 
ture  can  be  generated  by  synthesis  from  all  elements  of 
the  stiffness  matrix.  For  example,  consider  the  simple 
structure  in  Fig.  8  and  the  structural  element  in  Fig.  9. 
The  stiffness  matrix  for  the  element  is  [K]^ 

«2ct — nr~ 

k 


u 

FIG.  8  SIMPLE  STRUCTURE  MODEL 


FIG.  9  ELEMENT  NO.  I 

where  k..  denotes  the  force  at  coordinate  q.  when  a  unit 
load  is 'applied  at  coordinate  q, .  Elements3^  and  3  and 
the  corresponding  stiffness  matrices  are  shown  in  Fig.  10 
and  11, 
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FIG.  10  ELEMENT  NO.  2 
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FIG.  II  ELEMENT  NO.  3 

Then  the  generalized  stiffness  matrix  is: 


kll 

k12 

k13 

kl4 

0  0 

k21 

k22 

k23 

k24 

0  0 

k3i 

k32 

(k33  +  k33  +  k33) 

(k34 

+  k»4) 

k35  k36 

k4l 

k42 

(k44 

+  k44^ 

H5  He 

0 

0 

k54 

k55  k56 

0 

0 

k<i3 

k64 

k&5  k66 

0 

0 

k73 

0 

0  0 

In  most  situations  some  coordinates  are  required  to 
define  displacements  for  which  there  are  no  associated 
masses.  Obviously,  the  computational  acceleration  a  =  P/m 
will  become  infinite  if  m  vanishes.  To  avoid  this  unreason¬ 
able  result,  the  number  of  generalized  coordinates  must  be 
reduced  to  the  same  number  as  the  number  of  masses  in  the 
structure. 

First  it  is  necessary  to  distinguish  between  those 
displacements  in  {q}  associated  with  the  masses  and  those 
associated  with  zero  masses.  Thus, 


<q> 

where 


Ml 

l<4>  J 


{ q > *  =  displacements  for  which  masses  [M]*  exist 
{q>^  *  displacements  for  which  masses  are  zero 
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The  equations  of  motion,  Equation  (7),  are  now  rewritten 
using  the  above  notation. 


This  equation  is  equivalent  to  the  two  following  equations. 


+  [K]11{q}»  +  [K]12{q}°  =  {Q}# - (9) 

[K]21{q}*  +  [K]22{q}°  =  {Q}° - (10) 

The  second  of  these  equations  is  solved  for  {q}° 

{ q}°  =  CK]"|  ({Q}°  -  [K]21{q}») - - - (.11) 

Substituting  { q } u  into  Equation  (9)  gives  the  reduced  form 
+  CK]s{q}»  =  { P}  - (12) 

where 

[K]s  =  m±1  -  CK]12[K]“|[K]21 - (13) 

{F}  =  { Q}* -  [K]12[K]22{Q}° - (1*0 


Now  Equation  (12)  can  be  treated  to  get  the  displacements 
{ q > *  of  the  masses.  The  displacements  {q}*  are  then  sub-  Q 
stituted  into  Equation  (11)  to  obtain  the  displacements  { q > 
for  which  the  masses  are  zero. 

5 .  Numerical  Solution 

a.  Equations 

The  equations  of  motion  are : 

+  [K]s{q}»  =  {F}  - (15) 

This  is  a  system  of  second  order  differential  equations. 
These  equations  may  be  replaced  by  a  system  of  first  order 
equations.  As  a  simple  example,  the  second  order  equation 

=  f(t,y,^-  ) 
dt*  dt 

becomes  the  system 
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&  -  V 

dt 


^  =  f(t,y,v) 
dt 

The  two  functions  y(t)  and  v(t)  are  now  computed  simul¬ 
taneously.  Initial  conditions  such  as  y I  .  .  *  yn, 

dy  i  t  tA  u 

-*•  ■  vn  become  y 

dt  .  _  u 


0* 


vQ  so  that 


the  generalized  initial  value  problem  is  established. 


After  the  equations  of  motion  are  reduced  to  first 
order  differential  equations,  the  Runge-Kutta  procedure 
can  be  used  to  solve  them.  The  Runge-Kutta  method  is 
briefly  discussed  as  follows3: 

Let  the  equations  be 

V’  a  fi/t.yjv),  v*  *  f2(t,y,v) 

where 

y'  =  y  =  y(t),  V  «*  v(t),  V'  =  — 

dt  dt 

The  formulas 


t— 1 

= 

hfl(tn- 
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+ 

V 
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where  h  =  step  size 


may  be  shown  to  duplicate  the  Taylor  series  for  both  func¬ 
tions  up  through  terms  of  order  four.  For  more  than  two  si¬ 
multaneous  equations,  say  n,  the  extension  of  the  Runge- 
Kutta  method  parallels  the  above,  with  n  sets  of  formulas 
required  instead  of  two. 

b .  Digital  Computer  Code 


A  computer  program  code  has  been  developed  to  solve 
the  drop  impact  problem.  This  program  is  listed  in  the 
appendix.  Only  the  lumped  masses,  cushioning  forces, 
initial  displacement  and  velocity  and  member  proper¬ 
ties  such  as  flexural  rigidity  El,  spring  constant  k 
and  member  length  L  are  required  for  input  data.  The 
computer  program  automatically  generates  the  stiffness 
matrix  of  the  structure  and  rhen  solves  the  equations  of 
motion  by  the  Runge-Kutta  method,  The  displacement,  velo¬ 
city  and  acceleration  at  any  coordinate  at  any  time  may 
be  printed  out.  The  printed  out  displacement  is  relative 
to  the  position  where  impact  begins.  From  these  displace¬ 
ments  the  relative  displacement  of  any  two  points  can  be 
computed.  The  computer  code  is  sufficiently  flexible  to 
serve  as  a  design  tool. 

6.  Experimental  Program 

a.  Model 

To  evaluate  the  computational  procedure  a  small  model 
has  been  prepared  and  used  for  collecting  experimental  data. 
This  model  which  is  shown  in  Fig.  1  was  designed  so  it  would 
be  relatively  simple  to  represent  with  a  lumped  parameter 
system  and  at  the  same  time  to  have  at  least  10  degrees  of 
freedom.  The  schematic  drawing  of  the  model  shown  in  Fig.  2 
suggests  how  it  should  be  represented  by  lumped  parameters 
and  Fig.  3  shows  the  coordinates  and  the  degrees  of  freedom. 
As  indicated  previously  the  model  can  be  thought  of  as 
representing  an  M37  truck,  although  it  was  not  designed 
specifically  for  that  purpose.  Its  primary  purpose  is  to 
provide  a  specific  and  convenient  structure  for  use  in 
evaluating  the  procedure  developed  and  previously  described, 
for  computing  structural  response  to  impulsive  loading. 

The  cushioning  system  for  this  model  was  designed 
to  provide  zero  bending  moment  at  the  center  of  mass  of 
the  model.  All  cushioning  forces  are  directly  applied 
to  the  masses  except  forces  F0  and  Fc  which  are  applied^  - 
to  the  structure  tied  to  massB2.  This  loading  is  similar 


to  that  applied  to  the  M37  truck,  since  the  dead  load  in 
the  bed  of  the  truck  is  not  usually  cushioned  directly, 
but  is  cushioned  through  the  structural  members  that  support 
it.  Small  masses  5  and  6  are  attached  to  the  main  masses 
1  and  2  respectively  and  are  cushioned.  Mass  7  is  optional 
and  not  cushioned. 

Only  vertical  motion  is  considered.  This  is  an  over¬ 
simplification  of  the  real  situation  because  non-vertical 
motion  is  always  possible  due  to  wind  drift  and  system 
oscillation.  However,  since  the  vertical  motion  is  the  most 
important  factor  in  airdrop  cushioning,  it  is  essential  that 
it  be  analyzed  separately,  and  well  understood  before  at¬ 
tempting  a  solution  of  the  more  general  problem.  Lateral 
and  longitudinal  motions  of  the  model  are  minimized  by 
the  two  vertical  steel  restraining  rods,  which  can  be 
seen  in  Fig.  1.  Stability  of  the  model  is  also  improved 
by  using  clusters  of  3  parallel  springs  instead  of  one 
helical  spring  for  the  support  of  masses  3  and  4. 

t. .  Measurements 

Accelerations  are  measured  with  fluid  damped  resistance 
type  accelerometers  mounted  on  masses  1,  2,  and  3.  Dis¬ 
placements  are  measured  at  supports  A,  C,  and  D. 

To  measure  the  deflections  a  special  device*  is  used. 

This  device  which  is  shown  in  Fig.  12  consists  of  two  U-shaped 
thin  steel  springs  with  strain  gages  mounted  as  shown. 

The  legs  of  the  U-shaped  pieces  can  be  deflected  toward 
each  other  a  considerable  amount  without  the  development 
of  appreciable  force  or  stress,  and  the  strain  gage  output 
is  proportional  to  the  deflection.  Other  better  known 
transducer  types  could  not  be  used  because  even  the  small 
lateral  motion  that  sometimes  occurs  might  damage  the  transducer. 

c .  Parameter  Values 

Numerical  values  for  the  model  parameters  are 
M^^  -  20  lb 
M2  =  30  lb 
M3  *  =  10  lb 

Mp.  *  Mg  ■  0.6  lb 

»  2.3  lb 

*  designed  and  developed  by  Dr.  C.  H.  Yew,  The  University 
of  Texas  at  Austin. 
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FI6.  12  TYPICAL  DEFLECTION  CAGE 


K  =  600  lb/in. 

E  «  29  x  10^  psi  (steel  beam) 

I  -  0.3393  in.4 
Fa  =  247.0  lb 
FB  *  Fc  =  183.5  lb 
pD  -  Fe  ■  120  lb 

*  2.0  lb-in. -sec2 

p 

J2  3  2.0  lb-ln.-sec 
Paper  honeycomb  pad  sizes  are 
A  =  3  x  2^  x  3  Inches 
B  =  C  =  2^-x2^x3  Inches 
D=*E  =  2x2^x3  Inches 

The  pad  sizes  are  selected  to  provide  an  acceleration 
of  12g,  assuming  a  crushing  strength  of  4,000  psf  for 
the  honeycomb.  This  crushing  strength  was  determined 
experimentally  by  making  dynamic  loading  tests  on  samples 
having  approximately  the  dimensions  of  these  pads. 

7 .  Discussion  of  Results 

The  most  recent  measured  and  computed  results  are  shown 
in  Figs.  13  through  20.  In  comparing  these  results  it  should 
be  noted  that  there  is  a  difference  of  some  significance 
between  the  cushioning  forces  in  the  mathematical  model 
and  in  the  laboratory  model.  For  the  former  force  is  as¬ 
sumed  constant  from  the  moment  of  impact  until  it  vanishes 
at  the  instant  the  velocity  of  the  cushioned  mass  vanishes, 
and  from  that  time  on  it  remains  at  zero.  On  the  other  hand 
the  cushioning  force  in  the  laboratory  model  has  a  modest 
initial  spike  and  then  decreases  constantly  until  the  velo¬ 
city  of  the  cushioned  mass  vanishes.  At  that  time,  which 
incidentally  does  not  coincide  with  the  time  the  corres¬ 
ponding  velocity  in  the  mathematical  model  vanishes,  the 
force  does  not  drop  immediately  to  zero  because  the  cushion 
has  some  resiliency.  The  differences  in  the  crushing  forces 
are  reflected  by  differences  in  the  accelerations  shown  in 
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(b)  COORDINATE  q„ 


(c)  COORDINATE  q. 


FIG.  14  MEASURED  AND  COMPUTED  VALUES  OF  q,  ,  q0  AND  ql9 
WITH  3%  (2.3  1b.)  UNCUSHIONED  MASS  AT  q„ 

(DROP  HEIGHT  18  INCHES) 
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(0)  COORDINATE  q, 
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FIG.  20  MEASURED  AND  COMPUTED  VALUES  OF  q,  AND  qs 
WITHOUT  2.3  lb  (3%)  UNCUSHIONED  MASS  AT  q„ 
(DROP  HEIGHT  18  INCHES) 


Figs.  15  and  20.  The  computed  accelerations  for  shown 
in  Fig.  15 *  are  almost  identical  in  form  to  the  assumed 
cushioning  force,  while  the  measured  accelerations  closely 
resemble  the  general  form  of  the  crushing  force-time-re- 
lation  for  paper  honeycomb  usually  observed.  The  rather 
high  frequency  oscillations  which  appear  on  the  measure^ 
acceleration  record  are  also  somewhat  typical  of  the  oscil¬ 
lations  seen  on  dynamic  stress-strain  curves  for  paper 
honeycomb.  The  source  of  these  oscillations  can  not  always 
be  pinpointed.  In  some  cases  they  are  believed  to  typify 
the  way  honeycomb  actually  crushes,  in  others  they  are 
believed  to  represent  vibration  of  some  part  of  the  testing 
system.  In  these  records  the  frequency  is  much  too  high 
to  be  the  uncoupled  frequency  of  M,  and  its  supporting  spring, 
and.  too  low  to  be  associated  with  the  elastic  body  vibra¬ 
tion  of  the  mass.  .  Consequently  the  oscillations  are  attri¬ 
buted  to  the  crushing  characteristics  of  the  honeycomb 
cushions . 

The  oscillations  which  appear  in  the  displacement, 
and  both  the  measured  and  computed  accelerations  of  M,  all 
have  very  nearly  the  same  frequency,  approximately  200 
Hertz.  This  is  also  the  uncoupled  frequency  of  M-^  It  is 
somewhat  surprising  to  see  oscillations  of  a  mass  as  large 
as  M,  at  200  Hertz.  If  the  peak  accelerations  associated 
with^the  oscillations  of  the  mass  are  computed  assuming 
a  steady  state  harmonic  oscillation  they  are  found  to  be 
of  the  order  of  2000g  whereas  the  computed  accelerations  do 
not  exceed  peak  values  of  about  2g,  and  the  measured  values 
are  of  the  order  of  6g,  Thus  the  oscillations  on  the  dis¬ 
placement  record  are  not  believed  to  be  connected  with  ac¬ 
tual  motion  of  the  mass.  They  are  believed  to  be  a  result 
of  vibration  of  the  displacement  gage.  The  fact  that  they 
occur  at  nearly  the  same  frequency  as  the  oscillations  in 
the  acceleration  records  is  probably  Just  coincidence.  If 
the  oscillations  are  ignored,  and  the  actual  curve  replaced 
by  a  smoothed  curve  in  each  of  the  6  records  shown  In  Figs. 

13  and  14. the  measured  and  computed  displacements  are  seen 
to  be  in  quite  good  agreement  considering  the  differences 
in  the  crushing  forces  discussed  above.  Differences  between 
the  measured  and  computed  deflections  would  undoubtedly  be 
much  greater  if  it  were  not  for  the  relative  Insensitivity 
of  displacement  to  the  form  of  the  acceleration  curve. 

The  oscillations  which  appear  in  the  computed  acceleration 
record  for  shown  in  Fig.  20  have  a  frequency  which 

is  nearly  equal^to  the  uncoupled  frequency  of  Mp.  No  mass 
in  the  system  will  vibrate  at  its  uncoupled  frequency  since 
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there  Is  coupling  between  all  masses  in  the  system.  Con¬ 
sequently  the  observed  frequency  should  not  be  expected 
to  agree  with  the  uncoupled  frequency.  The  frequency  of 
the  oscillations  in  the  measured  acceleration  is  almost 
twice  the  frequency  seen  in  the  computed  accelerations. 

No  logical  explanation  for  this  discrepancy  can  be  of¬ 
fered  at  this  time. 

Although  the  computed  and  measured  values  of  the  dis¬ 
placements  and  accelerations  do  no*-  agree  precisely  they 
are  close  enough,  it  is  believed,  to  indicate  that  with  the 
same  crushing  force  in  the  two  models  the  results  could 
be  brought  into  close  agreement.  This  means  that  the  con¬ 
tinuous  laboratory  system  can  be  replaced  for  computations 
by  a  lumped  parameter  system  which  will  adequately  repre¬ 
sent  the  essential  features  of  the  motion  of  the  major 
parts  of  the  laboratory  model.  Hence,  a  vehicle  could  be 
represented  in  the  same  way. 

It  may  be  further  noted  that  the  addition  the 
small  mass  does  not  significantly  affect  the  displace¬ 
ments  and  accelerations  of  the  other  masses.  This  means 
that  only  the  major  masses  in  a  physical  system  need  to 
be  represented  in  the  lumped  parameter  model. 

8.  Conclusions 


Although  the  comparisons  between  computed  and  experi¬ 
mental  results  are  not  exhaustive,  and  there  is  a  signi¬ 
ficant  difference  in  crushing  forces  in  the  two  models,  the 
following  conclusions  are  believed  to  be  justified. 

a.  Displaceme  ;ts,  velocities,  and  accelerations  of 
the  principal  parts  of  a  laboratory  model  can  be  satis¬ 
factorily  predicted  using  a  lumped  parameter  matnematical 
model  and  a  numerical  computation  procedure. 

b.  If  predictions  can  be  made  for  a  laboratory  model, 
then  predictions  of  comparable  accuracy  are  possible  for 
actual  prototype  vehicles. 

c.  The  computed  program  is  simple  enough,  ard  the  com¬ 
putational  time  required  is  short  enough  to  make  the  pro¬ 
gram  a  practical  tool  for  design  purposes. 

d.  The  insignificant  effect  which  the  small  masses 
hu.ve  on  the  motions  of  the  larger  masses  in  the  system  in¬ 
dicates  that  models  of  prototype  vehicles  can  be  quite 
simple  and  still  give  an  accurate  indication  of  the  move¬ 
ments  of  the  different  parts  during  an  impact. 
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e.  Procedures  for  establishing  parameter  values  such 
as  mass,  and  stiffness,  for  lumped  parameter  models  to  re¬ 
present  actual  vehicles  are  not  indicated  in  the  results  of 
this  study.  These  procedures  remain  to  be  developed. 

9.  Recommendations  for  Further  Studies 

a.  Although  the  results  obtained  for  the  one  dimen¬ 
sional  model  were  not  entirely  conclusive,  further  work 
should  be  concentrated  on  a  two  dimensional  lumped  mass 
model.  For  the  theoretical  analysis  of  such  a  model  it 
is  necessary  only  to  modify  the  stiffness  matrix  as  a 
grid  structure.  A  grid  structure  would  not  only  be  a 
more  realistic  representation  of  a  real  vehicle  body,  it 
would  also  make  the  laboratory  model  more  stable.  This 
would  improve  the  accuracy  of  measurements. 

b.  The  damping  factor  which  was  not  considered  in 
this  report  should  be  investigated  in  further  studies. 

This  can  be  done  rather  easily  mathematically  but  the 
question  of  how  to  introduce  controlled  damping  into 
the  laboratory  model  requires  some  study. 

c.  An  actual  vehicle  for  which  measurement  data 
are  available,  or  can  be  obtained,  should  be  modeled 

so  the  procedure  can  be  tested  in  a  real  life  situation. 

d.  The  significance  of  the  parameters  which  are 
being  determined  should  be  carefully  reviewed  in  the 
context  of  vehicle  damage. 
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Flow  Chart  and  Computer  Program 
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FIG.  21  MAIN  PROGRAM  FLOW  CHART 
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NUMERICAL  SOLUTION  OF  A  LUMPED  PARAMETER  FOR  AIRDROPPED 
VEHICLES  BY  USING  RUNGE-KUTTA  METHOD. 

CDC  6600  t  FORTRAN  IV  ,  STORAGE  REQUIRED  30000  WORDS 
DIMENSION  SK ( 30,  30 )  .NOD ( 30 ) »U< 30 » 30 ) »C < 30 .30 ) »NV( 30 ) .  SQ<30) 
DIMENSION  DF(IO.IOO) .YJ(IOO)  .YJO(lOO) 

DIMENSION  YJ3 ( 30 ) »  SQ3i30) 

COMMON  N,S(30,30),GQA!30),WM(30)»WG<30) 

FORMAT (3110, F20. 2) 

FORMAT! 21 10,2 F 20. 4) 

FORMAT! 1615) 

FORMAT !8F10.3 ) 

FORMAT ( 16F5. 1 ) 

FORMAT  (  1H1  »*  BEAM  MBR  =  *  ,  I  5. 2X  » *SPR  ING  N0.  =  *  I5,2X,*TOTAL  JNT  N0.  =  * 
1  »I5»2X*#Es*,Fl2»2) 

FORMAT!I4,12F11.0/(4X,12F11.0)  ) 

INITIALLY  THE  PROGRAM  ASSIGNS  A  SET  OF  COORD.  ACCORDING  TO  THE 
ORDER  OF  JOINT  NUMBER  •  IT  BECOMES.  FOR  EXAMPLE,  Q ( 1 1 )  AND  Q ( 12 ) 
OR  THE  VERTICAL  DISPLACEMENT  AND  THE  ROTATION  AT  JOINT  NUMBER  6. 
READ  IN  NS  =NO.  OF  STRUCTURES 

Nll=NO.  OF  COORD.  WITH  ASSOCIATED  MASSES 

N22=NO.  OF  COORD.  WITHOUT  ASSOCIATED  MASSES 

NCH=  1  ,  IF  THE  ORDER  OF  COORD.  WILL  BE  REARRANGED 

NCH=  0  ,  NO  REARRANGEMENT  IN  COORDINATE  ORDER. 

READ!  5,5040)  NS 

RE AD ! 5 ,5040 )  Nll.N22.NCH 

JC=N1 1+N22 

READ  IN  MBR  =  NO.  OF  BEAM  MEMBERS 
NP  =  NO.  OF  SPRINGS 
JNT  =  TOTAL  JOINT  NUMBER 
E  =  MODULUS  OF  ELASTICITY 
READ ( 5  »  5010 )  MBR.NP.JNT.E 
WRITE  <6,6010)  MBR.NP.JNT.E 
G=386.04 
N=2*N11 

mt=mbr+np 

jm=jnt-np 

JC=2* JNT-NP  $  JNT 1  =  JC~ 1 
DO  5  1=1, JC 
DO  5  J=1,JC 
SK ( I , J ) =0. 

DO  100  1=1. MT 

READ  IN  MEMBER  PROPERTIES, 

JOINT  NUMBERS  CAN  BE  ASSIGNED  ARBITRARILY  IN  THE  ORDER  BUT 
SHOULD  BE  CONTINUOUS  FROM  1  THROUGH  THE  LAST  JOINT  NUMBER, 

NO  NUMBER  IN  THIS  INTERVAL  CAN  BE  OMITTED.  ALL  SPRING  JOINTS 
NOT  ATTACHED  TO  BEAMS  SHOULD  ALSO  BE  NUMBERED  FOLLOWING  THE 
BEAM  JOINT  NUMBERS. 

NJ1,  NJ2=  JOINT  NUMBERS  At  ELEMENT  ENDS 

A I =  MOMENT  OF  INERTIA  FOR  BEAM  MEMBER,  OR  SPRING  CONSTANT 
SL=  LENGTH  OF  BEAM  SEGMENT  , 

SL  =  0.0  FOR  SPRING 
READ ! 5,5020 )  NJ1 ,NJ2, Al »SL 
WRITE  <6,5020)  NJ1 »NJ2 » AI , SL 
IF  <  SL )  20,20,10 

GENERATE  STIFFNESS  MATRIX  OF  STRUCTURE 
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FA=E*AI/SL**3 

S1=12.*FA  $  S2  =  6.*SL*FA  $  S3=2 .*5L*SL*F  A  $  S4=2.*S3 

I F ( NJ 1-NJ2  )  11.11.12 

11  N 1  =  2*N Jl~  1  $  N3  =  2*NJ2-1 

GO  TO  16 

12  N1=2*NJ2-1  i  N3=2*NJ1-1 
16  N2  =  Nl  +  i  $  N4  =  N3  + 1 

SK(N1  ,N1)=SK(N1,N1  )+Sl  $  SK(N1,N2)=SK(N1»N2)-S2 

SK(Nl ,N3!=oK(Nl  .N3)-S1 

SK(N1.N4)“SK(N1,N4)-S2  $  SK ( N2 .N2 ) = SK ( N2 .N2 ) +S4 

SK(N2.N3)=SK<N2.N3)+S2 

SK ( N2  »N4 ) =SK ( N2  »N4 ) +S3 

SK(N3,N3)=SK(N3,N3)+S1  $  SK(N3»N4)=SK(N3»N4)+S2 

SK(N4»N4)=SK(N4,N4)+S4 

GO  TO  100 

20  IF(NJl-JM)  30.30,32 

30  N1=NJ1*2-1 

GO  TO  40 
32  N 1 =NJ 1+JM 

40  I F ( NJ  2- JM )  50,50.55 

50  N2=NJ2*2-1 

GO  TO  60 
55  N2=NJ2+JM 

60  CONTINUE 

SK ( N1 »N1 )=SK( N1 jNl J+AI  $  SK < N2 » N1 ) =SK ( N2 .Nl ) -A  I 
SK(N2,N2)=SK(N2.N2)+AI 
100  CONTINUE 
N1=JC-1 

IF(NCH)  280,280.210 

C  -  READ  IN  NOD ( I ) =NEW  ORDER  OF  COORD. 

C  -  ALL  INPUT  AND  OUTPUT  FROM  HERE  ON  ARE  REFERRED  TO  THIS 

C  -  NEW  COORDINATES  SYSTEM 

210  RE AD ( 5 , 5G40 )  ( NODI  I ) , I =1 , JC ) 

WRI TE ( 6 ,5040 )  (NODI  I ) , 1  =  1, JC ) 

!  =  JC 

220  J=NOD(I) 

SK ( I . 1 ) =SK( J. J) 

1  =  1-1 

I F ( I )  230,230,220 

230  DO  240  1=1, JC 

240  SK ( I , I ) =SK (1,1) 

DO  260  J= 1 ,N 1 
J1=J+1 
JP=NOD( J ) 

DO  260  I=J1,JC 
I P  =NOD ( I ) 

IF(IP-JP)  251,251.252 

251  SK( I . J)=SK( IP, JP) 

GO  TO  260 

252  SK( I.J)=SK( JP.IP) 

260  CONTINUE 

280  CONTINUE 

DO  270  J= 1 , N 1 
I 1=J+1 

DO  270  1=1 1, JC 
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SKI J.I )«SK( I.J) 

270  CONTINUE 

DO  140  I»1,JC 

140  WRITE (6,6020)  I  * ( SK( I » J ) » J* 1 » JC ) 

CALL  MIV(SK»U»N11»N22) 

DO  320  1*1  •  N 1 1 
DO  320  J=1*N22 

SUM=0. 

DO  310  Ka 1 »N22 
K1=K+N11 

310  SUM*SUM+SK ( I  *  K1 ) *U( K  » J ) 

C  < I ,J  )=SUM 
320  CONTINUE 

DO  340  1  =  1. Nil 

DO  340  J= 1 »N 1 1 
SUM=0. 

DO  330  K=1.N22 
K1=K+Nll 

330  SUM=SUM+C ( I»K)*SK(K1»J) 

340  S( I.J)=SK( I»J)-SUM 
DO  341  1  =  1, Nil 

341  WRITE  (6.6020)  I , ( S ( I »J ) » J=1 ,N1 1 ) 

339  CONTINUE 

C  -  READ  IN  H  =  TIME  STEP  SIZE  OF  INTEGRATION  IN  RUNGE-KUTTA  METHOD. 

C  H  =  0.0  OR  BLANK  FOR  STARTING  NEXT  STRUCTURE  OR 

C  LAST  CARD  OF  INPUT  DATA  DECK 

C  DT  =  TIME  INTERVAL  OF  PRINT  OUT 

C  TEND=END  OF  TIME  OF  CALCULATION 

READ ( 5 .5050 )  H.DT .TEND 
I F (H. EQ.  0.)  GO  TO  380 
READ  (5,5040)  LM.LW.LT.LV.LS 
NG0=0 

IF(LM)  343,343,342  * 

342  NG0=NG0+1 

READ( 5,5050 )  ( WM( I ) *  I  =  1 ,N1 1 ) 

C - READ  IN  W.4(I)=  MASSES  (IF  LM  GREATER  THAN  0) 

C  WG ( I ) =5TAT I C  REACTION  AT  CUSHIONED  COORD.  (IF  LW  .GT.  0) 

C  NT=TOTAL  NUMBER  OF  CUSHIONED  COORD.  (IF  LT  .GT.  0) 

C  NV( I )=CUSHIONED  COORD. 

C  V0=EQUI VALENT  FREE  FALL  VELOCITY  ( 1 F  LV  .GT.  0) 

C  YJO ( I ) =1 •  AT  ALL  VERTICAL  DISPLACEMENTS 

C  Y JO ( I ) =0.  AT  ROTATIONAL  COORD. 

C  SQ( I )=CUSHIONING  FORCES  (IF  LS  .GT.  0) 

343  WRI TE ( 6,5050)  ( WM ( I ) , I =1 ,N1 1 ) 

IF(LW)  345,345,344 

344  NG0=NG0+1 

READ (  5,5050)  (WG( I ) . 1=1 ,JC ) 

3<k>  WR  t  TE  ( 6 ,5050  )  ( WG(  I  )  » I  =1 ,  JC  ) 

IF(LT)  347,347.346 

346  NG0=NG0+1 

READ( 5,5040)  (NT.  (NV ( I  ) ,  I =1 .NT > ) 

347  WRI TE( 6,5040 )  (NT.  (NV ( I ) , I  =  1  .NT ) ) 

IF(LV)  349,349,348 

348  NG0=NG0+1 

READ( 5,5050 )  VO 
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WR I TE ( 6  *5050 )  VO 

READ (  5,5060)  ( Y J0( I ) , I  =  1 , JC  ) 

DO  351  I*1,JC 
Y  JO ( I ) =VO*Y  JO (  I  ) 

351  Y J 3 ( I >  =V JO (  I  ) 

GO  TO  353 

349  DO  352  I=1»JC 

352  Y  JO (  I ) =  YJ3 (  I  ) 

353  WR I TE ( 6*5050 )  (YJO< I ) ♦ 1=1, JC) 

IF(LS)  355,355,354 

354  NG0=NG0+1 

350  READ  (5,5050)  ( SQ ( I ) ♦ 1 = 1 ♦ JC  ) 

DO  359  I*1,JC 

359  SQ3(I>*5Q(I> 

GO  TO  357 

355  DO  356  1=1, JC 

356  SQ(I) =SQ3 ( I  ) 

357  WR  I  T fc  (  6  ,  *•050  )  (  SQ  (  I  )  ,  I  =  1 ,  JC  ) 

IF(NGO)  380,380,358 

358  CONTINUE 

DO  21  I«1*JC 

J=  I  *2 

YJ( J-1)=0. 

21  YJ(J)=YJO(I) 

DO  370  1  =  1, Nil 
SUM=0. 

DO  360  J= 1 »N22 
Jl= J+Nl 1 

360  SUM=SUM+C( I , J)*SQ( J1 ) 

GQA ( I )=SQ( I ) -SUM 

370  CONTINUE 

WRITE! 6,5050 )  ( GQA ( I )  ,  I  =  1 ,N1 1 ) 

PRINT  385 

385  FORMAT!  1K1 , // »4X ,*COORD*, 10X»*DI SP. * , 1 5X ,*VELOCI TY* , 12X ,*ACC . 
CALL  RKAMSB! 3,1. 0E-7, 1.0,1. OE-1 0,1. 0E-8, 0.0) 

ALPHA=0.0 

K5=2*JC 

DO  141  1=1, K5 

141  Y JO (  I ) =Y J I  I ) 

390  OMEGA=ALPHA+DT 

IF (OMEGA. GT. TEND)  GO  TO  339 
150  CALL  RKAMSB ( 1 , ALPHA ♦ OMEGA *H,S5E»YJ) 

PRINT  14, OMEGA 

14  FORMAT ! 20X  »3H  T=,  F15.6) 

CALL  DERFCN  (  OMEGA , YJ , 1 ,DF > 

DO  15  1=2, N, 2 

11=1/2 

12=1-1 

15  WRI TE ( 6,6001  )  1 1 , Y J ( I  2 ) » YJ ( I ) ,DF ( 1 , I  ) 

N I =N1 1+1 

DO  410  I  — N I  ,N 
SUM=0 • 

DO  400  J= 1  ,N 1 1 
J 1= J+2-1 

400  SUM  =  SUM+SK ( I »  J ) *Y  J ( J1 ) 
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GQA( I )=SQ( I)-SUM 
410  CONTINUE 

DO  430  1  =  1, N22 
SUM=0. 

DO  420  J=1,N22 
J1=J+N11 

420  SUM=SUM+U( I , J)*GQA( J1  ) 

I 1=N1 1+1 
K1  =  I 1*2-1 
K2=K1+1 
YJ(K1)=SUM 

Y  J ( K2 )  =  ( SUM- Y JO ( K 1 ) J/DT 
DF(1,I1)=(YJ(K2)-YJ0(K2))/DT 
WR I  TE ( 6 ,600 1 )  II, SUM  ,YJ(K2) ,DF< 1, I  1) 

430  ONTINUE 
MQ=0 

DO  515  1  =  1, NT 
J=NV ( I ) 

K  =  J*2 

460  IF(YJO(K)*YJ( K)>  470,515*515 

470  I F  C  SC ( J ) )  480,490,490 

480  SQ( J ) =0.0 
GO  TO  500 
490  SQ(J)=-WG(J) 

500  MQ= 1 
515  CONTINUE 

IF (MQ  )  550,550,520 

520  DO  540  1  =  1, Nil 
SUM=0 • 

DO  530  J=1 ,N22 
J 1=J+N11 

SUM=SUM+C( I,J)*SQ(J) 

530  CONTINUE 

540  GQA( I )=SQ< I J-SUM 

550  CONTINUE 

ALPHA=OMEGA 
GO  TO  390 
380  CONTINUE 
NS=NS-1 

IF (NS )  26,26,1 

6001  FORMAT (5X» I4,4E20,5) 

26  STOP 
END 
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SUBROUT 1NE  RKAMSB ( MODE . A1 .A2 . A3. A4, Y J ) 

DIMENSION  X( 10) ,Yi 10.100) »YJ{ 100) »0{ 10.100 ).DF( 10.100) »XS(  10)  . 
1  W( 10) .B( 10.100) .A (10) *BP( 10).BC( 10) .PHI ( 100) .YS( 10*100) 
DOUBLE  Y.PHI.YS 

COMMON  N.S<  30.30) .GOA ( 30) »WM< 30) »WG( 30) 

GO  TO  (2*2*1) *MODE 

1  CONTINUE 

C  -  COEFFICIENTS 

W(  l)=WI4)=l./6.  S  W(2)=W(3)=l./3.  $  KK  =  4 

A(2)=A(3)=*5  S  A  (  4  )  =  1 . 

B(2»1)=B(3»2)=.5  5  B ( 3 . 1 ) =B ( 4 , 1 ) =B ( 4.2 ) =0.  $  B(4.3)=l. 

IQ  =  4  $  FCT  =  19./270. 

BP ( 1 ) =-9 • / 24.  S  BP(2)=37./24.  S  BP ( 3 ) =-59. /24.  $  BP ( 4 ) =55 . /24. 
BC ( 1 ) = 1 • / 24.  *  BC<2)=-5./24.  S  BC ( 3 ) =19. /24.  $  BC(4)=9./24. 

C -  TOLERA  ces 

HMI N  =  A 1  S  HMAX  =  A2  S  EM  I N  =  A3  S  EMAX=A4  S  ISKIP  =1 

I F  (  YJ(l).EQ.O.)  ISICIP  =  0 

RETURN 

2  CONTINUE 

c  -  RKAMSUB  ENTRY  POINT 

CALL  DERFCN(Al.YJ.l.DF) 

ALPHA  =  A1  S  OMEGA  =  A2  $  H  =  A3 
I F  (  MODE.EQ.  2  )  10  •-  KK 

I QM1  =  10-1  $  I0P1  =  10+1  $  I STP=0  S  S I GN= 1 . 

I F (  H.LT.O.  )  SIGN  =  -1. 

X (  1  >  =  ALPHA 
DO  3  1=1, N 

3  Y  (  1  *  I  )  =  Y  J  (  I  ) 

4  MM  =  1  $  I FLG  “  0 

5  KCGUNT  =  0 

6  M  =  MM  S  MM  =  M  +  1 
I F (  MM.GT.IQPJ  )  MM  =  1 
X ( MM )  =  X ( M )  +  H 

TEST  =  OME ) A  -  X(MM)  $  TEST1=TEST/0MEGA 
IF  QUOTIENT  OVERFLOW  7,8 

7  TEST1  =  TEST 

8  I F (  ABS(TEST1).LT.1.0E-10  >  GO  TO  12 

I F (  SIGN*TE3T1  )  9,12.13 

9  T EST 2  =  OMEGA  -  X(M) 

I F (  MODE.EQ. 2  )  GO  TO  11 

I F (  I  SKIP. EQ. 0 .OR.S IGN* TEST 2.LT.HMIN  )  GO  TO  99 
H  =  TEST2/IO 

I F (  S IGN*H.LT *HMI N  )  GO  TO  11 
M  =  M-l 

I F (  M.EO.O  )  M= I  OP  1 

X  (  1 )  =  X  (  M  ) 

DO  10  1=1, N 

10  Y(  1.1)  =  Y ( M  ,  I  ) 

GO  TO  4 

11  H  =  TEST2  $  IFLG  =  0 
X ( MM )  =  X(M>  +  H 

12  ISTP  =  1 

13  XJ  =  X ( M) 

DO  14  1  =  1 v  N 

1^*  Y  J  (  I  )  =  Y  (  M,  *  ) 
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I F (  MODE.EQ. l.AND.M.EQ.IQ.OR.IFLG.EQ.l  )  GO  TO  32 

C  -  RUNGE-KUTTA  PROCEDURE 

DO  25  K=1»KK 
I F <  K.EQ.l  )  GO  TO  22 
XJ  =  X ( M )  +  H*A(K) 

DO  15  I  =  1»N 
15  PHI ( I  )  =  0. 

KM3  =  <  -  1 
DO  20  1=1, N 
DO  19  J=1,KM1 

19  PH  I  ( I  )  =  PHI (  I  )  +  H*B ( K  » J )*D ( J, I  ) 

20  YJ(  I  )  =  Y  <  M, I )  +  PHI (  I  ) 

22  CALL  DERFCN(XJ»YJ»K»D) 

I F (  IFLG.EQ.2.0R.K.NE.1  )  GO  TO  25 
DO  23  1=1, N 

23  DF ( M»  I  )  =  D( 1  *  I ! 

25  CONTINUE 

DO  27  1=1, N 
27  PHI ( I)  =  0. 

DO  30  1=1, N 
DO  29  K= 1 »  KK 

29  PH  I  ( I  )  =  PHI (  I  )  +  H*W ( K) *D ( K » I ) 

30  Y ( MM , I )  =  Y ( M , I  )  +  PHI (  I) 

I F (  ISTP.EQ.l  )  GO  TO  100 
GO  TO  6 

c  -  ADAMS  PREDICTOR-CORRECTOR  PROCEDURE 

32  CALL  DERFCN' X J, YJ»M»DF ) 

DO  33  1=1, N 

33  PHKII  =  0. 

DO  34  K=1 » IQ 

J  =  K  +  KCOUNT 

I F <  J.GT.IQP1  )  J  =  J  -  I QP1 
DO  34  1=1, N 

34  PH  I  ( I )  =  PHI ( I )  +  H*BP(KJ*DF( J,I ) 

DO  35  1=1, N 

35  Y J ( I )  =  YIM.I )  +  PHI ( I ) 

XJ  =  X ( MM ) 

CALL  DERFCN(XJ,YJ,MM,DF) 

DO  43  1=1, N 

43  PHI(I)  =  0. 

DO  44  K=1,IQ 

J  =  K  +  KCOUNT  +1 
I F (  J.GT.IQP1  )  J  =  J  -  IQP1 
DO  44  I =1 , N 

44  PHI1I)  =  PHI (  I  )  +  H*BC(K)*DF( J,  I  ) 

DO  45  1=1, N 

45  Y ( MM, I )  =  Y ( M , I )  +  PHI (I) 

C  - —  SINGLE-STEP  ERROR 

DLTMX  =  0. 

DO  54  1=1, N 

DLT  =  ABS(  Y ( MM , I )  -  YJ(I)  ) 

I F  <  DLT. LE. DLTMX  )  GO  TO  54 
DLTMX  =  DLT  S  IDLT  =  I 
54  CONTINUE 

TEST  =  DLTMX /Y ( MM, IDLT ) 
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IF  QUOTIENT  OVERFLOW  55*56 

55  SSE  =  ABS (  FCT*DLTMX  )  S  GO  TO  60 

56  SSE  =  ABS (  FCT*TEST  ) 

60  CONTINUE 

C  -  ERROR  ANALYSIS 

I F  {  ISMP.FQ.O  )  GO  TO  90 

I F {  EMIN. LT. SSE. AND. SSE. LT.EMAX  )  GO  TO  90 
I F (  SIGN*H.GT.HMIN.AND.SIGN*H.LT.HMAX  )  65*89 

65  I F  <  SSE.GT.EMAX  )  66*70 

66  H  =  H/2. 

I F (  SIGN*H.LT.HMIN  )  89*68 

08  I F (  IFLG.EQ.O  )  GO  TO  4 
M  =  M-l 

I F (  M.EQ.O  )  M  =  I QP 1 
X( 1)  =  X l M ) 

DO  69  1=1, N 

69  Y  (  1. 1  )  =  YIM,  I  ) 

GO  TO  4 

70  I F (  ISTP.FQ.A  )  GO  TO  100 
73  H  =  2.*H 

I F (  SIGN*!  XIMMl+H-OMEGA  )  )  77*75,75 

75  H  =  H/2.  S  GO  TO  90 

77  I F (  SIGN*H.GT.HMAX  )  89,78 

78  IBK  =  IQ/2  +  1  S  L  =  0 
DO  82  K  =  1 » I BK 

J  =  I BK  -  K  +  1  $  M  =  MM  -  L 

I F {  M.LE.O  )  79,80 

79  M  =  M  +  I QP 1  $  L  =  0  S  MM  =  M 

80  XS(J)  =  X ( M i 
DO  81  1=1, N 

D( J, 1  )  =  DF(M,  I) 

81  Y  5 1  J  ,  I  >  -  Y ( M , I) 

82  L  =  L  +  2 

DO  85  K= 1 » I BK 
X ( K )  =  XV  K) 

DO  85  1=1, N 
DF ( K , I )  =  D ( K  » I  ) 

35  Y (K,  I  )  =  YS ( K  *  I  ) 

MM  =  IBK  $  IF  LG  =  2  5  GO  TO  5 

89  PRINT  511  $  GO  TO  100 

90  I F (  ISTP.EG.l  )  GO  TO  100 

IFLG  -  1  $  KCOL'NT  =  KC0UNT  +  1 

I F (  KCOUNT.EQ. IQP1  )  KCOUNT  =  0 
GO  TO  6 
99  MM  =  M 
100  CONTINUE 

c -  RKAMSUB  EXIT  POINT 

A2  =  X ( MM )  $  A3  =  H  %  A4  =  SSE 

DO  105  1=1, N 
105  Y J ( I  )  =  Y ( MM , I  ) 

5x1  FORMAT! 1H1,5X,25H  STEP  SIZE  OUT  OF  BOUNDS  ) 

500  RETURN  $  END 
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SUBROUTINE  DERFCN  <XJ»YJ,M»DF) 

COMMON  N,S(30,30),GQA(30)»WM(30)»WG(30) 
DIMENSION  DF( 10,100) ,YJ( 100) 

N l=N/2 

DO  10  I  =  1 .N . 2 
10  DFIM.I )=YJ( 1+1) 

DO  30  I*2»N»2 
SUM=0.  S  11=1/2 
DO  20  J= 1 »N1 
J1=J*2-1 
J2=J1+1 

20  SUM=SUM+S< I1»J)*YJ< Jl> 

30  DF ( M» I )  =  ( GQA( I  1 ) -SUM ) /WM( 1 1 ) 

RETURN  $  END 


SUBROUTINE  MIV(SK.U.Nll.NM) 

DIMENSION  A( 30.30) »U( 30.30) .SK( 30.30) 
DO  9001  1  =  1, NM 
I  1= I+Nll 
DO  9001  J= 1 »NM 
J 1=J+N1 1 

A( I . J)=SK( I1.J1) 

U ( I » J  )  =0. 

IF  (I.EQ.J)  U(I.J)=1.0 
9001  CONTINUE 

EPS=0. 0000001 
DO  9C15  1=1, NM 
K  =  I 

IF  (I-NM)  9021.9007,9021 
9021  IF  ( A ( I  ,  I ) -EPS  )  9005,9006,9007 

9005  I F ( -A ( I » I ) -EPS )  9006,9006,9007 

9006  K=K+1 

DO  9023  1 , NM 

U( I ,J )=U( I , J)+U(K,J) 

9023  A( I, J)=A(I »J)+A(K,J> 

GO  TO  9021 

9007  DI V=A { I , I ) 

DO  9009  J= 1  ,NM 
U( I » J)=U( I ,J) /DIV 

9009  A(I»J)=A(I ,J)/DIV 
DO  9015  MM=1 »NM 
DELT=A(MM,I) 

IF  ( ABS(DELT)-EPS)  9015,9015,9016 
9016  IF  ( MM- I )  9010,9015,9010 

9010  DO  9011  J= 1 » NM 

U ( MM , J ’ *U ( MM , J ) -U ( I ,J)*DELT 

9011  A { MM, J ) =m( MM, J ) -A( I,J)*DELT 
9015  CONTINUE 

RETURN  $  END 
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•  The  accelerations  and  displacements  in  a  complex  structure  subjected  to 
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independent  coor Jinates  through  which  the  motion  of  these  lumped  masses  are 
uniquely  described,  and  by  writing  equations  of  motion  in  terms  of  these 
coordinates,  a  set  of  equations  is  derived  which  represents  the  motion  of 
any  part  of  the  model  during  impact.  Using  the  Runge-Kutta  numerical  method, 
and  a  digital  computer,  these  equations  are  solved. 

A  physical  model  of  '.he  lumped  parameter  system  was  built  and  cusioned 
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